The random Blume-Capel model on cubic lattice: first order inverse freezing in a 3D 

spin-glass system 
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We present a numerical study ol the Blume-Capel model with quenched disorder in 3D. The phase 
diagram is characterized by spin-glass/paramagnet phase transitions of both first and second order 
in the thermodynamic sense. Numerical simulations are performed using the Exchange-Monte Carlo 
algorithm, providing clear evidence for inverse freezing. The main features at criticality and in the 
phase coexistence region are investigated. The whole inverse freezing transition appears to be first 
order. The second order transition appears to be in the same universality class of the Edwards- 
Anderson model. The nature of the spin-glass phase is analyzed by means of the finite size scaling 
behavior of the overlap distribution functions and the four-spins real-space correlation functions. 
Evidence for a replica symmetry breaking-like organization of states is provided. 

PACS numbers: 



I. INTRODUCTION 

The so-called inverse transition (IT) is a reversible 
transformation occurring between phases whose entropic 
contents and symmetries are in the inverse order relation 
relatively to standard transitions. The case - already hy- 
pothesized by Tammann more than a century agoi - of 
"ordering in disorder" taking place in a crystal solid that 
liquefies on cooling, is generally termed inverse melting. 
The IT phenomenon also includes the transformation in- 
volving amorphous solid phases - rather than crystal - 
as that of a liquid vitrifying upon heating. In this case 
the term inverse freezing is somewhat used in the liter- 
ature: both phases are disordered but the fluid appears 
to be the one with least entropic content. The reason for 
these counter intuitive phenomena is that a phase usually 
present only at high temperature happens to exist also 
in peculiar patterns such that its entropy is actually less 
than the one of the phase normally considered the most 
ordered one. 

Inverse transitions in their most generic meaning (i.e., 
both thermodynamic or occurring by means of kinetic 
arrest) have been detected in the last years in a number 
of different materials and between phases of rather differ- 
ent nature. The first example was the transition between 
liquid and crystal phases of helium isotopes He 3 and He 4 
at low temperature i 2 - A more complex and recent exam- 
ple is the polymer poly(4-methylpentene-l) - P4MP1 - in 
which a crystal polymer melts as the temperature is de- 
creased, or the pressure increased. By means of exhaus- 
tive measurements by Differential Scanning Calorime- 
try (DSC) and X-ray diffraction the phase diagram of 
P4MP1 has been experimentally determined by the group 
of RastogipT— showing evidence for both an equilib- 
rium inverse melting, between a crystal phase (tetrag- 
onal or hexagonal, depending on the pressure) and a 
fluid phase, and a non-equilibrium IT between the hexag- 
onal crystal and a glassy phase. Another extensively 



studied instance is a molecular solution in water, com- 
posed by a-cyclodextrine (evCD) and 4-methylpyridine 
(4MP) mixed in given molecular ratios, investigated by 
means of neutron scattering, X-ray diffraction, DSC and 
rheometric measurements The "solid" is in this case 
a sol-gel porous system formed by an ordered network 
of molecules of aCD-water-4MP filled with liquid 4MP, 
melting down decreasing temperature at constant aCD 
concentration. Eventually, another important polymeric 
example is methyl-cellulose solution in water, undergoing 
a reversible inverse sol-gel transition! 16 i 17 For such sys- 
tem, a careful analysis of the behavior of the microscopic 
components across the transition has been performed in 
literature^ and, therefore, it turns out to be particularly 
important for the modelization proposed in the present 
work, as we will see in the following. 

Apart from polymeric and macromolecular substances, 
in the last years ITs showed up in many other different 
contexts. Inverse melting from an ordered lattice to a dis- 
ordered vortex phase takes place, e.g., for the magnetic 
flux lines in a high temperature superconductor^ A gas 
of atoms at zero temperature passes from superfluid to 
insulator as the lattice potential depth is increased* 2 ^ 
Furthermore, in the framework of nanosystems, the re- 
versible transition of an isotropic liquid into an ordered 
cubic phase upon heating has been detected experimen- 
tally in ferromagnetic systems of gold nanoparticles i 2 1 1 22 

In this work we stick to a definition of IT as the 
one put forward by Tammann) 1 - a reversible transition 
in temperature at fixed pressure - or generally speaking, 
at a fixed parameter externally tuning the interaction 
strength, such as concentration, chemical potential or 
magnetic field - from a solid high temperature phase to an 
isotropic fluid (or a paramagnet, for magnetic systems) 
low temperature phase. Generalizing to non-equilibrium 
systems one might address as IT also those cases in which 
the isotropic fluid is dynamically arrested into a glassy 
state. This occurs, e.g., for the crystal-glass transition 
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in the cited P4MP1 as pressure is not too large 3 or in 
molecular dynamics simulations and mode-coupling com- 
putations of attractive colloidal glasses . 23 > 24 

In this definition IT is not an exact synonym of reen- 
trance. Indeed, though a reentrance in the transition line 
is a common feature in ITs, this is not always present, 
as, e.g., in the case of aCD^ 7 -^ or methyl-cellulose 1 ^ so- 
lutions for which no high temperature fluid phase has 
been detected. Moreover, not all re-entrances can be 
seen as signatures of an IT to a completely disordered 
isotropic phase. In liquid crystals, ultra-thin films and 
other materials, phases with different kind of symmetry 
can, actually, be found that are separated by reentrant 
isobaric transition lines in temperature - cf., e.g., Refs. 
[25U29I ] - but no melting occurs strictu sensu. Also re- 
entrances between dynamically arrested states, aperiodic 
structures or amorphous solids of qualitatively similar 
nature, like liquid-liquid pair a 30 ' 31 are not considered as 
IT, since an a-priori order relationship between the en- 
tropic content of the two phases is not established and 
it cannot be claimed what is inverse and what is "stan- 
dard". For the same reasons also re-entrances between 
spin- glass (supposed at equilibrium, that is, considered 
as a thermodynamic phase) and ferromagnetic phases - 
as, e.g., in Refs. 32.33] - hardly fall into the IT cate- 
gory. Eventually, re-entrances in parameters other than 
temperature are also not taken into account as inverse 
melting/freezing transitions. 

A thorough explanation of the fundamental mecha- 
nisms leading to the IT would need of a microscopic anal- 
ysis of the single components behavior and their mutual 
interactions as temperature changes across the critical 
point. Due to the complexity of the structure of poly- 
meric chains and macromolecules involved in such trans- 
formations, a clear-cut picture of the state of the single 
components is often not available. For the case of the 
above mentioned methyl-cellulose, Haque and Morris 18 
proposed that chains exist in solution as folded bun- 
dles in which hydrophobic methyl groups are packed. As 
the temperature is raised, the bundles unfold, exposing 
methyl groups to water molecules and, thus, causing a 
large increase in volume and the formation of hydropho- 
bic links eventually leading to a gel. The polymers in 
the folded state are thus inactive (or far less interact- 
ing than those in the unfolded state) but also yield a 
smaller entropic contribution than the unfolded ones. As 
the chains start to unfold because of thermal noise they 
change to an interacting state thus enforcing bonds with 
other chains and condensing in a gel. 

Theoretical modeling for IT is starting to develop but 
is still on its first, often uncorrelated, steps and con- 
sists, at the better, in heuristic reproductions of the 
phenomenon^ - — Looking, in particular, at the transi- 
tion between an amorphous 'frozen' phase and a fluid 
(i.e., paramagnet), recently spin-glass models with spin-1 
variables have turned out to effectively represent systems 
in which the transformation is driven by entropic effects. 
In these cases inverse freezing has been studied in the 



mean-field approximation! °' 

We also mention that with the help of this class of mod- 
els, the connection between entropy driven phase reen- 
trance and shear thickening can also be tackled 4 ^ and, 
furthermore, a generalization of the spin-1 variable to a 
composition of "fast" and "slow" variables- 4 coupled to 
two different thermal baths allows for studying anoma- 
lous latent heat in out of equilibrium transitions^ 3 - 

In the present work we will consider the Blume-Capel 
(BC) model^ 4 - with quenched disorder: a spin-glass model 
on a 3D cubic lattice with bosonic spin— 1 variables 
(sj = — 1,0, +1). Under the assumption that the inter- 
play between inactive and interactive states of a micro- 
scopic component is at the ground of the eventual IT, 
bosonic spins can approximate the folded/unfolded con- 
formation, S = representin g th e inactive state, S = ±1 
the interactive one, cf. Ref. 38] for a more comprehen- 
sive discussion. We will focus on the random version of 
the BC model introduced by Ghatak and Sherrington 45 
(GS) in order to study the effects of the crystal-field in 
a spin glass - e.g., (Ti 1 - x V x ) 2 03 displays anisotropic 
spin glass behavior in function of x 4^ The mean-field 
solution in the Full Replica Symmetry Breaking (RSB) 
schem o 40 ' 47 " — predicts a phase diagram with a second 
order transition line between spin-glass (SG) and para- 
magnetic (PM) phase ending in a tricritical point where 
a first order phase transition line starts and a phase co- 
existence region appears^ Furthermore, the first order 
transition is characterized by the phenomenon of IT i 40 i 49 
the low temperature phase is PM with a lower entropy 
than the SG phase and the transition line develops a 
reentrance. 

In the original (ferromagnetic) BC model a 44 i 50 how- 
ever, no IT was observed in the mean-field approxima- 
tion, nor in finite dimension studies^ 1 - - — and in presence 
of quenched disorder a recent study on a 3D hierarchi- 
cal lattice by means of renormalization group theory in 
position spaced 5 - provides no evidence for a low tempera- 
ture tricritical point or a PM/SG reentrance, contrarily 
to what is predicted by mean-field theory. 

Moreover, we found in the literature only one finite 
dimensional system with quenched disorder undergoing 
a standard first order phase transition in finite dimen- 
sion: the 4-Potts glass studied by Fernandez et al. in 
Ref. [56[ . In that work a first claim has been made that 
first order phase transition exists in 3D systems also in 
presence of quenched disorder, though the randomness 
tends to strongly smoothen the transition into a second 
order one. This transition is driven by the temperature 
and by the degree of dilution of the Potts glass bonds. 
Though in numerical simulations changing, e.g., the pres- 
sure, the bond dilution, or even the relative probabilities 
of the random bond values, cf., e.g., Ref. [57j, is techni- 
cally equivalent, the latter are complicated to control in 
a real experiment and require the preparation of several 
samples with different microscopic properties. The study 
of a conceptually simpler model, satisfactorily approach- 
able with standard simulation techniques, might help in 
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validating the assessment of the existence of first order 
phase transitions in random systems. 

Motivated by the above considerations we have, thus, 
studied the existence of inverse freezing in the 3D disor- 
dered BC model with nearest-neighbor interactions and 
the nature of the "frozen" (or, rather, "blocked") phase. 
We present hereafter the results of our investigation by 
means of Monte Carlo numerical simulations. Some re- 
sults about the critical behavior have already appeared 
in a recent letter^ 

In the present manuscript we first introduce the model, 
in Sec. II, and in Sees. Ill and IV we define the numer- 
ical techniques employed to study continuous and dis- 
continuous phase transitions in finite size (FS) systems. 
In Sec. V we recall the Exchange Monte Carlo method, 
else called Parallel Tempering (PT)^^ In Sec. VI, we 
present our results about the phase diagram of the model 
and its critical behavior both along the continuous tran- 
sition and in the coexistence region related to the first 
order transition. The main features of the organization 
of states in the SG phase in finite dimension (i.e., below 
the upper critical dimension for our model) is studied 
in Sec. VII, where we perform a systematic study of 
the properties of the overlap distribution functions and 
of the four-spins correlation functions in space. Finally, 
Sec. VIII reports our conclusions. 



significant order parameter that drives the transition is 
the density p of magnetically active (|sj| = 1) sites: 

i 

The apex J recalls us that the values of the parameters 
depend on the particular realization of disorder ({Jij})- 
Useful information about the equilibrium properties of 
the system can be obtained from the knowledge of the 
following probability distribution functions (pdf) 



P(q s ) ee Pj(q s ) = (d (q„ - qi J) )) (6) 

P( qi ) ee Pj{qT) = (S (qi~ql J) )) (7) 
P(p) EE pJ(p) = (6(p-pV))) (8) 



where . . . denotes the average over quenched disor- 
der and (...) the thermal average. Though the den- 
sity probability distribution is known to be self-averaging 
(limjv->-oo Pj,n{p) = Pj(p)), this does not hold for the 
overlap distributions Pj{q s ,i)r^ f° r which 86 



P(q*,i) = Pj(q s .i) # lim Pj,n(Q:i) ( 9 ) 



II. MODEL AND ORDER PARAMETERS 

We consider the following Hamiltonian 

where (ij) indicates ordered couples of nearest-neighbor 
sites, and s$ = — 1,0, +1 are spin— 1 variables lying on 
a cubic lattice of size N = L 3 with Periodic Boundary 
Condition (PBC). Random couplings are independent 
identically distributed as 

P(Jij) = \8(Jij-l) + \s(JiJ + l) (2) 

The field D is usually called crystal-field and it plays the 
role of a chemical potential for the empty sites s = 0. 
We will, therefore, refer to D invariably as chemical po- 
tential or crystal field in the following. We simulate two 
real replicas of the system and define their site and link 
overlaps, i.e., the order parameters usually characterizing 
the SG transition, as 

«P - ^E-PM 95 (3) 

i 

where 3 is the dimension of the space. If a thermody- 
namic phase transition occurs, with latent heat, the most 



Through the study of the pdfs, as function of the ex- 
ternal thermodynamic parameters, we can identify the 
PM /SG transition and discriminate between first and 
second order phase transitions. As an instance, if a first 
order phase transition takes place, the density pdf P(p) 
displays a double peak due to the coexistence of the PM 
and SG phases. Moreover, by means of the overlap pdfs 
we can investigate the nature of the SG phase. 

III. FINITE SIZE SCALING FOR CONTINUOUS 
TRANSITIONS 

In order to infer the details of the critical behavior 
from numerical simulations of finite size systems, a fun- 
damental quantity (in zero external magnetic field) is 

C ^ = jf E i^s+r) 2 (10) 
s 

with r — (r x ,T y ,r z ). In terms of space-dependent over- 
laps, q r = Sr, C4 can be written as 

Ci{r) EE ^E^p9p+r)l2 (11) 
P 

- 1 V(s (1) s (1) )i(s (2) s {2) )o 

— Z_^\ b P b p+r/l\ b P *p+r/2 

V 

where (. . .)x2 stays equivalently for the thermal average 
((• ■ -)i)2 or ((• • -)2)i over the two replicas independently. 
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This is the four-spins correlation function, and the infor- 
mation it provides can be exploited in different ways to 
identify the existence of a second order phase transition 
for finite size systems and probe the thermodynamic be- 
havior in the SG phase. 

A conventional way to identify a second order phase 
transition is to look at the behavior of a correlation 
length-like scaling function defined as 



_ / drr 2 C 4 {r) <91ogC* 4 (fc) 



JdrC 4 (r) 



dk 2 



(12) 



fc 2 =0 



where 



C 4 (fe) 



(2tt) 3 



d :i re- lkr C 4 (r) 



On a 3D cubic lattice, the above defined correlation 
length becomes^ 



1 



( CM 

4 sin 2 \ \ (7 4 (fei; 



- 1 



(13) 



where k\ = \k\\, k\ = (2%/L,0, 0) is the minimum wave- 
vector of the lattice and = (0,0,0). In the thermody- 
namic limit, a second order transition is characterized by 
a diverging correlation length, at critical temperature T c , 
whose Finite Size Scaling (FSS) behavior is the same as 
in Eq. iPjl^Zl 

Another relevant observable is the SG susceptibility 



Xsg = lV) = i 3 c 4 (o) 



(14) 



diverging at the PM/SG transition as L — > oo. Because 
of FS, though, £ c and xsg cannot diverge in numerical 
simulations, although inside the critical region a remark- 
able property of the critical phenomena survives: scale 
invariance. Indeed, we can define a FS "critical" temper- 
ature Tg, function of the size L, as the temperature at 
which the above mentioned observables do not (or only 
slightly) depend on the size: 

I = &(|) =t(L 1/v (T-T c L )) (15) 
Xsg^- 2 = X (|) = X(L 1/U (T - T C L )) (16) 

The values T£ at which £ C /L at different L cross each 
other are the FS respective of the critical temperature. 
The latter can. thus, be estimated by FSS in the L — > oo 
limit. 

A further size independent observable at criticality is 
the so-called Binder parameter: 



9(L,T) 



1 



94 
(92)' 



(17) 



with q n = ((qi' J,> ) n ). It measures the deviation of P(q) 
from a Gaussian distribution as the SG phase is ap- 
proached. Since q 4 and q 2 scale with L in the same way, 
g does not depend on L at T c . 



A. Quotient method 

Denoting by 0(T, L) a generic observable diverging at 
critical temperature T c as L —> oo, and considering two 
sizes L, L' whose scale ratio is s — L'/L, we can look at 
the scaling of quotient 



0(T, sL) 
0(T,L) 



s)+0(Z-«,L-») (18) 



where Fo is a universal FSS function and uj is the power 
of the subleading FS corrections. Through the scaling 
Ansatz, (1181) one, thus, introduces a class of universal 
functions Fo that are size-independent in the critical re- 
gion. Given two observables O and R displaying scale 
invariance, this allows for plotting Fo vs Fr for several 
values of L: if the data collapse in a universal scaling 
function, the scaling Ansatz is verified and FSS methods 
can be trusted to evaluate the critical exponents. We 
will analyze in the present manuscript the behavior of 

F t> f xsg and F g- 

In order to estimate the critical exponents we can, 
thus, use the so-called quotient method^ based on the 
observation that at = T* , the correlation lengths in 
systems of sizes L and sL (in L units) are equal: 



sUT* c ,L) = 

Indent For an observable O scaling as t x ° (t 
in the thermodynamic limit, it holds: 



(19) 



T/T c - 1) 



Q(t;,sl) 

0(T*,L) 



(20) 



where the dependence through in the correction term 
is neglected because, in the critical region, £ c ^> L. For 
a SG we can obtain the exponents v and r\ by means of 
the FSS of the quotients of <9^£ and xsg, scaling, respec- 
tively, with exponents 



i 



(2-77)1/ 



These relations hold if the Ansatz (IT51) is verified,— i.e. 



if Fo is a size-independent scaling function for several 
values of L and sL. 



IV. CHARACTERIZATION OF A FIRST 
ORDER TRANSITION 

Besides the second order transition, the random BC 
mean-field model also shows a tricritical point beyond 
which a true first order phase transition with non-zero 
latent heat occurs and from which a region of coexistence 
of PM and SG phase departs4S The slope of a first order 
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line is given by the Clausius-Clapeyron equation that, for 
our model reads 



dD _ s PM - ssg _ As 
dT p PM - Psg Ap 



(21) 



where D plays the role of the pressure. The equilibrium 
transition line changes slope in a point where the entropy 
of the two coexistence phase is equal As = (Kauzmann 
locus^i). 

In order to tackle the identification of a first order 
phase transition in a 3D finite size system from numerical 
simulation data we sketch in the following four methods 
to estimate critical (and spinodal) points. 
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A. "Equal weight" estimate 

The transition takes place at the point at which the 
configurations belonging to the SG phase and those be- 
longing to the PM phase have the same statistical weight, 
i.e., they yield identical contribution to the partition 
function of the single pure phase. Else said, the free 
energies of the two coexisting phases are equal. In this 
statistical mechanical framework, the FS transition line 
D C (L,T) can be obtained as the locus of points where 
the two phases are equiprobable, i.e., the areas of the 
two peaks are equals 

rpo /•! 

/ dpP(p) = / dpP(p) (22) 

Jo J pa 

where po € [ppm '■ Psg] such that P(po) = (or mini- 
mal next to the tricritical point). A way to numerically 
determine the transition point is, thus, to compare the 
areas under the peaks of the distributions, cf. Eq. ([5]). 



B. Maxwell "Equal distance" estimate 



FIG. 1: Graphical sketch of the equal distance (dashed/green 
arrows) and equal area methods (dotted/blue areas) for FS 
systems undergoing first order phase transition for a given 
instance (L = 6, T = 0.5). 

in the (p, D) plane at a given temperature and deter- 
mine the value of D c as the one whose corresponding p c d 
value along the D(p) FS curve is equally distant from 
both ppm (-D c ) and psg{D c ), cf. Fig. Q] 

Pcd{D c ) = ~ [p PM {D c ) + p SG (D c )} (23) 

C. Maxwell "Equal area" estimate 

Alternatively we can determine D c as the value at 
which, cf. Fig. [Tj 

/>D C rDpM rDpM 

/ PPM (D)dD+ PsG (D)dD= p(D)dD 

Jd sg JD c J D SG 

where Dqq and -Dpm are arbitrary, provided that they 
pertain to the relative pure phases. 



There are other two methods to determine a first order 
transition in finite systems, based on the Maxwell con- 
struction. If we are in the coexistence region the curve 
D(p) for the system of size L will display a sort of plateau 
around some D = D* ~ D^: in a very small interval of 
D the density changes very rapidly. In the case of a pure 
state, instead, the D(p) curve has a far smoother behav- 
ior. In order to estimate the critical point we need to 
extrapolate the behavior of D(p) for the pure phases in- 
side the region of coexistence. To this aim we perform 
two fits exclusively based on the points outside the spin- 
odal points: for D < D^ G for the SG phase (Dsg(p)) 
and D > D s p p M for the PM phase (D PM (/o)). We will call 
Ppm,sg(-D) the inverse of the curves -Dpm.sg(p) extrap- 
olated from the data points pertaining to the pure PM 
and SG phases, respectively. The curve p(D) will denote 
the inverse of D(p). 

In this way we can make a Maxwell-like construction 



D. Symmetric distribution estimate 

Defining the skewness of the density probability distri- 
bution as 



C((p» 



((p-(p)) 3 ) 

((p-(p))2)3/2 



(25) 



in Ref. [66] the critical point was estimated as the point 
at which the double peaked distribution is symmetric. 
Since the skewness of P{p) can be precisely computed 
this estimate does not suffer, e.g., of the arbitrariness of 
the choice of po. In the thermodynamic limit, indeed, in 
the phase coexistence region both peaks of P(p) should 
be Dirac delta distributions and equal weight would be 
equivalent to a symmetric bimodal distribution. We will 
show the outcome of this analysis in our system for dif- 
ferent cases and compare it with the previous ones. 



6 



V. EXCHANGE MONTE CARLO ALGORITHM 
IN T AND D 

Wc have simulated our spin-1 model using the parallel 
tempering (PT) algorithm, replicating several copies of 
the system both at different temperatures and at different 
values of the external field D. For the PT in temperature, 
the swap probability of two copies at temperature T and 
T + AT is: 

P S wa P (A/3) - min [1, exp{A/3AH}] , (26) 

with A/3 = 1/(T + AT) - 1/T; whereas the swap proba- 
bility in the chemical potential reads 

P Bvrap (AD) = min [1, exp{f3ADAp}} (27) 

We used the latter implementation in trying to identify 
the reentrance of the transition line in the T, D plane. 
Since, however, the transition turns out to be first or- 
der in the whole region of inverse freezing, the PT al- 
gorithm must be handled with caution. Indeed, at the 
transition Ap is discontinuous implying the vanishing of 
-Pswap(AD) around the critical point for a given fix probe 
AD. In order to overcome this problem we have used a 
varying AD, smaller in the candidate coexistence region 
and larger and larger outside. An instance of this kind 
choice is represented in Fig. [2] 

For very large sizes, though, this would require a too 
precise a-priori knowledge of the transition lines and the 
method could not be applied with a reasonable success. 
On the other hand, the FSS effects appear to be almost 
nonexistent at the first order phase transition, so that 
probes at larger sizes were actually not necessary. In 
Tabs. HI |TI]we report our simulation parameters for PT 



in temperature and in crystal field, respectively. 
Thcrmalization has been checked in three ways. 

1 . We have verified the symmetry with respect to zero 
of the site overlap distribution function Pj(q s ) for 
single random samples (cf, Fig. [3]). In absence of 
an external magnetic field this must be symmetric 
for any choice of {Jij} realization. 

2. We have looked at the i-log behavior of the energy 
and we have considered as thermalized those sys- 
tems in which at least the last two points coincide 
within the error, cf. Fig. [3] This means that at 
least the last half of the data in MCS can be used 
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TABLE I: Simulation parameters of the parallel tempering 
in temperature: number of samples 2000, Monte Carlo Steps 
(MCS), number of thermal bath N T spaced by AT =,0.02 or 
0.025. 
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FIG. 2: Values of the chemical potential D for the replicas 
in the PT simulation exchanging systems at different D. The 
parameters refer to the simulated L = 12 system at T = 0.2. 
The dashed (green) line is the estimate of the FS critical value 
D C (L, T) estimated by means of the equal weight method (cf. 
Sec. llVfl . Inset: Chemical potential intervals AD vs. D in 
log scale for the same instance. 





L 


6 


8 


10 


12 


15 


0.2 


Din 


1.99 


1.999 


2.00392 


1.981 


1.981 




AD in 


0.002 


0.0006 


0.00027 


0.003 


0.003 




N D 


21 


21 


37 


37 


37 




MCS 


2 15 


2 17 


2 18 


2 20 


2 20 


0.3 


Din 


2.0034 


2.026 


2.0212 


2.0256 


2.028 




AD in 


0.002 


0.001 


0.00037 


0.003 


0.00025 




N D 


21 


21 


21 


31 


31 




MCS 


2 15 


2 17 


2 17 


2 17 


2 18 


0.4 


Din 


2.05 


2.06 


2.057 


2.06 


2.062 




AD in 


0.003 


0.002 


0.0007 


0.00085 


0.0006 




N D 


21 


21 


21 


31 


31 




MCS 


2 15 


2 17 


2 17 


2 17 


2 18 


0.5 


Din 


2.06 


2.06 


2.06 


2.026 


2.026 




AD in 


0.01 


0.01 


0.01 


0.008 


0.008 




N D 


21 


21 


21 


37 


37 




MCS 


2 15 


2 17 


2 17 


2 18 


2 is 



TABLE II: Simulation parameters of the parallel tempering 
in D. Number of disordered samples: 1000 
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FIG. 3: Instances of thermalization checks. Top: Pj(q s ) and 
Pj(-Qs) for an arbitrary sample at D = 0, T = 0.5, L — 16. 
Bottom: average energy versus time (in MCS) in log scale, 
t = log(#MCS)/log2. 

for computing statistical ensemble averages. 

3. we have checked the lack of variation on logarithmic 
time- windows of all considered observables (e.g., £ c 
and xsg) on at least two log points. 



VI. NUMERICAL RESULTS. 

A. Second order phase transition and universality 

In Fig. @]we present the T-behavior of £, C /L for differ- 
ent values of D = 0, 1, 1.75, 2 and 2.11. In the first four 
cases the curves at different L clearly cross, yielding evi- 
dence for a non zero T C L . From a FSS T C L = T c °° + aL- b 
we can, thus, extrapolate the critical temperature in the 
thermodynamic limit. The T C L crossing points between 
L — 2L curves are reported on column 3 of Tab. IIIIl The 
FSS estimates are reported in columns 2 and 5 of Tab. 
IIV1 where L/L' couples are chosen both as L/2L (col. 2) 
and as contiguous in the series L = 6, 8, 10, 12, 16, 20, 24 



(col. 5). Analogous, though less precise estimates, can 
be obtained by studying the behavior of the Binder cu- 
mulant g, cf. Eq. (fl~7|) . Applying both the quotient and 
the conventional FSS methods we can, eventually, obtain 
two estimates for the critical exponents. 

Before applying these methods, though, we must check 
if we can exclude cross-over effects as the chemical po- 
tential D is varied due to FS. Since we are in presence of 
a tricritical point, signaled, among others, by the weird 
behavior of £ C /L a t D = 2.11, cf. Fig. HI we should 
control how it influences the results as it is approached 
along the continuous transition line increasing D. In the 
mean-field approximation, indeed, at the tricritical point 
the coefficient of the fourth order term in the SG free 
energy action goes to zero and the sixth order term be- 
comes relevant for the critical behavior, as shown in Ref. 
[47l | . This is a typical behavior of Blume-Emery-Griffiths- 
Capel-like system a 67 ' 68 that might hinder the determina- 
tion of the critical behavior in the neighborhood of the 
tricritical point for sizes that are "not large enough" . 
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6 
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1.02(4) 


1.35(1) 


2.31(1) 


5.1(1) 
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8 


- 16 


0.99(6) 


1.31(2) 


2.58(5) 


5.1(1) 


-0.36(3) 




10 
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1.0(1) 


1.35(4) 


2.3(1) 


5.2(1) 


-0.39(4) 




12 


- 24 


0.98(9) 


1.33(2) 


2.43(7) 


5.1(1) 


-0.35(4) 




00 


1.01(1) 




2.34(3) 




-0.36(1) 
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L 


— L' 


Tc(s) 


Qded s ) 


u(s) 


Q x (s) 


77(B) 




6 


- 12 


0.894(9) 


1.32(1) 


2.51(4) 


4.9(3) 


-0.29(9) 


1.0 


8 


- 16 


0.895(9) 


1.396(6) 


2.08(1) 


4.8(4) 


-0.26(1) 




10 


-20 


0.877(9) 


1.271(7) 


2.89(2) 


5.1(5) 


-0.33(1) 




12 


- 24 


0.86(1) 


1.35(1) 


2.29(2) 


5.1(5) 


-0.3(1) 




00 


0.88(1) 




2.45(1) 




-0.31(2) 
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L 


— L' 


Tc(s) 


Qd d s ) 


u(s) 


Qx(s) 






6 


- 12 


0.715(7) 


1.41(1) 




4.7(5) 




1.75 


8 


- 16 


0.679(9) 


1.37(1) 


2.12(4) 


4.8(5) 


-0.3(1) 




10 


- 20 


0.67(1) 


1.34(3) 


2.4(1) 


5.0(6) 


-0.3(1) 




12 


- 24 


0.68(1) 


1.38(1) 


2.11(3) 


4.9(5) 


-0.3(1) 




00 


0.69(1) 




2.20(3) 




-0.30(1) 


D 


L 


— L' 


Tc(s) 


Qd d s ) 


u(s) 


Q x (s) 


r)(s) 




6 


- 12 


0.593(7) 


1.59(4) 




9.5(9) 


* 


2.0 


8 


- 16 


0.569(8) 


1.47(3) 


* 


18(2) 


* 




10 


-20 


0.54(1) 


1.37(3) 


* 


16(2) 


* 




12 


- 24 


0.54(1) 


1.34(4) 




10(2) 


* 




00 


0.56(1) 











TABLE III: Critical temperature and exponents are calcu- 
lated with QM: for D = 0.00, D = 1.00 and D = 1.75, 
through a FSS analysis of the values of Qa 3 j(s) and Q x { s ) 
for s = L'/L = 2. Cells with * mean that quotients are com- 
puted on sizes too small to significantly represent the asymp- 
totic behavior with L. 
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—I 0.9 1 1.1 1.2 0.8 0.9 1 0.6 0.7 0.8 
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T 

FIG. 4: Scaling functions £ C /L vs. T for different values D. For D = 0, 1, 1.75, 2 (L = 6, 8, 10, 12, 16, 20, 24) evidence for a 
continuous phase transition is found in the region of scale invariance. At D = 2.11 [L = 6, 8, 10, 12) no crossing is observed 
and, at low T, £ C /L — > 0. 

I 



To estimate and control FS effects we use the scaling 
methods introduced in Sec. IIIII and compare different 
universal FSS functions. In Fig. [5] we plot the Binder 
parameter g vs. the rescaled correlation length £ C /L at 
all simulated values of the chemical potential D both for 
a small (L — 6, top) and a large (L — 20, bottom) sys- 
tem. In the top plot one can easily observe that as the 
tricritical value of D is approached (2.05 < D 3c < 2.11) 
for L = 6 the curves do not overlap with each other sig- 
naling an apparent lack of universality. At large enough 
sizes, instead, all curves are superimposed (bottom plot 
of Fig. [5j L — 20), demonstrating that universality holds 
along the whole continuous transition line and that, be- 
cause of strong FS effects, a crossover occurs and the 
analysis limited to (or including also) too small sizes can 
hinder the prediction of the asymptotic behavior. 

The same effect is clearly shown in Fig. [6] where the 
size dependence of spin-glass susceptibility at criticality 
is shown. As D increases towards D^ c there appears to 
be a crossover in the scaling moving from small to large 
sizes and induces wrong asymptotic values of the critical 
indices. We, thus, did not make use of the small val- 
ues of L for D ~ D^c, namely L — 6 at D = 1.75 and 
D = 6,8, 10 and 12 at D = 2, to interpolate the values of 



the critical exponents, as they induce a wrong estimate 
as the limit L — > oo is performed. 

As a test for the eye, in Fig. [71 we display g vs. £ c /£ 
for all D and L values employed for our FSS analysis. 
Without the smallest sizes near the tricritical point, uni- 
versality appears quite tidy. In Fig. [8] we parametrically 
plot the universal FSS functions F%, F XSG and F g , cf. Sec. 
IIIII vs. £ C /L, as well, for the same simulated systems. 

The critical values of the temperature and the expo- 



D 


T c 


V 


n 


T c 


V 


V 


0.00 


1.01(1) 


2.34(3) 


-0.36(1) 


1.0(1) 


2.5(2) 


-0.37(2) 


1.00 


0.88(1) 


2.45(1) 


-0.31(2) 


0.8(1) 


2.6(5) 


-0.31(2) 


1.75 


0.68(2) 


2.20(3)* 


-0.30(1)* 


0.6(1) 


2.6(6) 


-0.30(4) 


2.00 


0.56(1) 


t 


t 


0.5(1) 


2.3(2) 


-0.31(2) 



TABLE IV: Critical temperature and exponents calcu- 
lated via QM Qafj^(s,T c (s)) and Qxsg( s >^ c ( s )) ( c °1 s - 2,3 
and 4) and via standard FSS analysis of the behavior of 
logd£c(L,T c (L))/d/3 and log X sg(L, T c (L j) (cols. 5, 6 and 7). 
*: estimated through QM without L = 6. f: not estimated 
by QM. 
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FIG. 5: Universal scaling function g vs. £ c /£ at L = 6 (left) 
and L = 20 (right) for all simulated D values. At small size 
the curves do not fall on top of each other as D is too near 
the tricritical point D = 1.75, 2, whereas at large size their 
critical behavior appears to be the universal for all D. 
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D=0 ,2-11=2.37(1) 
D=1 ,2-11=2.31(2) 
D=1.75, 2-11=2.30(4) 
D=2 ,2-11=2.31(2) 
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FIG. 6: xsg vs. L in log-log plot for D = 0, 1, 1.75, 2. For D = 
2.00, near the tricritical point, a cross-over is evident from 
small (L < 12) to large (L > 12) sizes. The quotient method 
does not yield reliable estimates because of a crossover in the 
scaling functions in the range of probed sizes. 




0.4 0.5 
l(L,T)/L 

FIG. 7: Binder parameter g(L,T c (L)) vs. £ C (L,T C (L))/L 
in the critical region for different D and L. Values of L = 
6, 8, 10, 12 for D = 2 and L = 6 for D = 1.75 are omitted. 
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FIG. 8: Scaling behavior of F XSG , F% and F g (top to bottom) 
vs. £ C (L,T)/L, cf., respectively, Eqs. CO}, (Q3]), and 
dHJ), at D = 0,1, L/2L = 6/12,8/16,10/20,12/24 and at 
D = 1.75, L/2L = 8/16, 10/20, 12/24. 



Model 


V 


V 


SG 3D bd G9 


2.22(15) 


-0.349(18) 


SG 3D bd™ 


2.53(8) 


-0.384(9) 


EA 3D2- 


2.15(15) 


-0.337(15) 


EA 3D^ 


2.00(15) 


-0.36(6) 



nents r] and v are shown in Tab. HVI both for the QM and 



TABLE V: Critical Indices of EA models in literature 



the canonical FSS methods. Due to the FS cross-over no 
interpolation was possible with QM at D — 2. We, thus, 
provide only one estimate for the indeces. 

As one can see, comparing with estimates of critical 
exponents summarized in Tab. [V] the system appears 
to be in the same universality class of the Edwards- 
Anderson model (corresponding to the D = — oo limit 
of our model) 2L22r2& 

In Fig. Iwe also plot £/L at D = 2.05 and D = 2.11 
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for L = 6, 8, 10, 12. In the first case we obtain a T c — 
0.553(7), though no analysis of the critical exponents can 
be performed because of FS effects. In the latter case no 
evidence is found for a second order phase transition, cf. 
Fig. |4j As T < 0.5 is approached, moreover, £ even ap- 
pears to scale weaker than L. We will see in the following 
why this comes about. 



B. First order phase transition 

Across a second order transition the system undergoes 
a transformation from a PM pure phase to a SG pure 
phase. As far as the density distribution P(p) is con- 
cerned, a pure phase corresponds to a single-peaked dis- 
tribution. As two peaks appear, the system exists both 
in PM (low p) and SG (high p) coexisting phases and 
we are in the neighborhood of a first order phase transi- 
tion. In FS systems the peaks are not delta-shaped but 
become sharper and sharper as L increases. At finite L, 
thus, P{p) is a good order parameter that drives the first 
order kind of transition: varying D, T, the system un- 
dergoes a transition with a discontinuous jump in p and 
the "thermodynamic" average values of p are obtained 
by looking at the peaks of its distribution. 

In Fig. [S] we show the behavior of the density dis- 
tribution through the first order transition at T = 0.4. 
The FS first order transition points can be determined 
with the four methods mentioned in Sec. IIV[ as we will 
show below. The spinodal lines at given L are estimated 
by looking at the D values at which a secondary peak 
arises. Since the region of phase coexistence corresponds 
to an inverse freezing transition, we performed PT sim- 
ulations at finite T, changing D. Indeed, in our model, 
we will see that the first order transition line displays 
a reentranc o i due to the existence of a "fluid" (PM) 
phase with an entropy lower than the one of the glassy 
phase. 

For what concerns the estimate of D C (T) the method 
of equal weight introduced in Sec. IIV1 cf. Eq. (|22[) works 
quite well for data collected at T < 0.4, because the two 
peaks are very well separated as soon as they appear, 
cf. Fig. [9l and the estimate is robust against reason- 
able changes of po (see inset of Fig. [§]). As T increases 
towards the tricritical value, however, the PM and SG 
values of the density approach each other. At T = 0.5, 
cf. Fig. [TTJJ we thus have the problem that the distribu- 
tions of the densities of the two phases are overlapping 
also for the largest simulated size. In that case, seen the 
arbitrariness of choosing po, we actually determine the 
transition point as the D value at which the peaks have 
the same height. This is a rough estimate but yields no 
difference w.r.t., e.g., fitting the two peaks separately and 
computing the areas under the interpolating curves. In 
Tab. I VII we report for all simulated sizes and tempera- 
tures the estimated values of the critical points obtained 
by this method, together with the spinodal points. 

These results can be cross-checked using the methods 
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FIG. 9: Density distribution P(p), L = 15, across the coex- 
istence region at T — 0.4: two peaks develop at ppm and 
pSG- As D increases the thermodynamically relevant phase 
(lowest free energy) passes from SG to PM in a first order 
phase transition. The dominant phase corresponds to the one 
with larger probability, i.e., larger integral of the peak. As 
the peak at psG vanishes the system is in a purely PM phase. 
Inset: Pis(p) on t/-Log scale. 
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FIG. 10: P(p) in the coexistence region at T — 0.5 and L = 
15. The two peaks are hard to distinguish and the coexistence 
region is rather narrow. 

based on the Maxwell construction, cf. Sec. IIVI and Fig. 
[TTT The pure phase behaviors Dpm,sg(p) are interpo- 
lated in the coexistence region by a polynomial fit on 
those points for which no double peak is present in the 
P(p). At any given L we look at the value D = D c such 
that (equal distance) 

PsgP>c - p(P>c) = p(D c ) - ppm(D c ) 

and at the value of D = D c at which the areas between 
D c and D(p) to the left and to the right of their crossing 
point are equal, i.e., 

AA(D C ) - / dp'(D(p r ) -D c )- dp'(D(p r ) - D c ). 

JO Jp{Dc) 

(28) 
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T 


D c 


n PM 
^SP 


n SG 
^SP 




T 


Dc[P{p)\ 


D c [pcd] 


D C [AA=Q)} 


D c [( = 0] 


0.2 


2.0031(1) 


1.9833(2) 


2.024(1) 
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2.0031(1) 


2.0033(2) 


2.0031(2) 


1.991(2) 
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2.032(3) 


2.015(1) 


2.043(5) 
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2.032(3) 


2.031(2) 


2.030(1) 


2.020(2) 
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2.060(1) 


2.046(2) 


2.092(5) 
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2.060(1) 


2.060(1) 


2.058(1) 


X 


0.5 


2.106(1) 


2.097(4) 


2.143(4) 
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2.106(1) 


2.103(3) 


2.102(1) 


X 



TABLE VI: Results of the first order phase transition: a fine TABLE VII: Evaluation of the first order critical point with 
tuning of the parameters {A} is needed in order to establish the method of equal weight (col. 2), equal distance (col. 3), 
the critical values D c , D SP and D SG - equal area (coL 3) and zero skewness - 
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FIG. 11: Maxwell construction at T = 0.2 in the p, D plane. 
L = 6,8, 10, 12, 15 from bottom to top as p > 0.6 (right). The 
almost vertical lines at the small and large density sides are 
the interpolated pure phase (PM left, SG right) behaviors. 
As an instance the critical D values for the equal distance 
construction are plotted at L — 6 (lower horizontal line) and 
at L — 15 (higher horizontal line). 
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FIG. 12: Estimates of D c vs. 1/L for L = 6, 8, 10, 12, 15 at 
T = 0.2, 0.3, 0.4, 0.5 (bottom up) obtained by the four meth- 
ods described in Sec. IIVI 



is zero. 

We, eventually, compute the skewness of double peaked 
P{p) as D changes, looking at the D = D c point for 
which Q(p{D c )) — 0. Since the two peaks of P(p) at 
finite size appear to be of different shape (SG broader, 



PM narrower), cf. Figs. [9j [TOJ the point at which the 
skewness is zero appears to be slightly different from the 
D c values computed with the previous three methods. In 
Fig. [T^] we plot at different temperatures the FS values of 
D C (L) with the four methods. The equal weight methid 
and the two Maxwell construction methods yield consis- 
tent results. For T = 0.4, 0.5 the estimate of D c by the 
symmetric distribution method displays a growing behav- 
ior in 1/T that does not allow for a consistent L — > oo, 
cf. Fig. QT] first and second panel from top, whereas 
at lower temperature, where the interpolated thermody- 
namic limit is stable the value is smaller than the other 
estimates. 

Summarizing, in Tab. IVIII we report the estimates of 
the first order critical point obtained by means of the 
four methods. 



C. Phase diagrams and inverse freezing 

Phase diagrams are plotted in Fig. [13] In the D, T 
plane we observe a pure SG phase at low T and D < 2. 
Increasing the temperature the continuous transition to 
the pure PM phase is denoted by a full line connecting 
the five numerical estimates of T c obtained by simula- 
tions at D = 0,1,1.75,2 and D = 2.05. We found no 
evidence for a continuous phase transition at D — 2.11. 
Beyond (D,T) = (2.05,0.53(2)) a tricritical point is 
placed. Beside changing to a first order transition, for 
lower T also a reentrance in the T C (D) line occurs. The 
warmest first order point for which we have an esti- 
mate is (D,T) = (2.109(2), 0.5). In Fig. ia detail 
of the phase coexistence region is plotted (inside the 
grey-dotted lines). In the inset of Fig. [13] we plot the 
(p,T) diagram. Below T — 0.53(2) no pure phase exists 
with an average p in between the dashed-grey curves. 

The inverse freezing takes place between a SG of 
high density to an almost empty PM (e.g., at T = 0.4, 
in the coexistence region D G [2.046(2) : 2.092(5)], 
Psg — 0.52 and ppm — 0.03). The few active sites 
do not interact with each other but only with inactive 
neighbors and this induces zero magnetization and over- 
lap. The corresponding PM phase at high T has, instead, 
higher density (e.g., p PM (D = 2,T = 0.6) = 0.4157(2), 
Ppm{D = 2.11, T = 0.6) = 0.596(2)) and the param- 
agnetic behavior is brought about by the lack of both 
magnetic order (zero magnetization) and blocked spin 
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case also the spinodal lines are reported (dashed). Inset: T, p 
phase diagram. 
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FIG. 14: Detail of inverse freezing region, interpolation of 
transition line D c (oo,T) (dotted), spinodal lines (dashed). 
The error bars are the FSS of the minimal interval in T and D 
at each L needed to identify the crossings in £ C /L curves (for 
continuous transitions) or compare the areas under Pn (p) for 
first order transition. In the top inset we show the latent heat 
|As|/T along the first order transition line. 



configurations (zero overlap). 

Using Eq. (|2~T|) . from the knowledge of Ap and the 
numerical estimate of dD / dT we are able to evaluate the 
latent heat employed in the transition, that we plot as a 
function of temperature in the top inset of Fig. Q3] 



VII. NATURE OF THE SG PHASE. 

The SG phase of the disordered BEG model, in mean- 
field regime, shows the same features of the Sherrington- 
Kirkpatrick model: in order to obtain a stable thermo- 
dynamics the Full RSB scheme is neededi 47 ' 48 On the 
other hand, out of the limit of validity of the mean-field 



regime, it is still unclear if the properties of SG phase are 
in agreement with the RSB scenario. The low T, D phase 
is characterized by a pure spin-glass phase and what this 
phase consists of in terms of statistical mechanic states 
is the subject of the following analysis. Three cases are 
contemplated in the literature. 

Droplet theory: it exists only one SG state (plus 
its symmetric spin-reversed) and, therefore, the over- 
laps between states in different replicas cannot fluctuate 
among different disordered samples and the distributions 
are delta-shaped* 7 - 7 - The four-spins correlation function 
in position space r = {x,y,z) should tend to a plateau 
C4 ( I r I ) = <?ea> f° r l ar g e enough |r|, that becomes longer 
as T decreases towards T c . 

Trivial- N on- Trivial (TNT) scenario: equilibrium 
states are many and non-trivially organized (i.e., q s fluc- 
tuates from sample to sample), but the excited states 
are droplet-like (i.e., the qi overlap, sensitive to inter- 
faces, fluctuates less and less as the size grows). This im- 
plies that P(q s ) is broad and non-trivial, whereas P(qi) is 
delta-shaped^ 7 - 8 - Since excitations are trivial, the expected 
behavior of 64(2;, y, z) is the same as for the droplet the- 
ory. 

Replica Symmetry Breaking (RSB) theory: many 
states characterize the SG phase, with space-filling exci- 
tations; both distributions are, thus, broad, with a com- 
plex structure i 61 ' 79 The correlation C±(x, y, z) is expected 
to decay continuously to zero (the minimum squared 
overlap for the present system, in absence of an exter- 
nal magnetic field) at all T i 81 ' 82 

First we will consider the overlap distribution func- 
tions, cf. Eqs. (©-([7]), since, in the spin glass phase 
(T < T c ), the site and the link overlap distributions - 
P(q s ) and P{qi) - can be used as hallmarks to discrimi- 
nate among different theories for finite dimensional spin 
glasses. In the next section we will analyze the four spin 
correlation functions. 

In order to see whether P(q s ) is trivial or not we need 
to estimate if, for growing sizes its support does shrink 
to a unique value, the Edwards- Anderson parameter <7ea 
or it remains finite. In our case, in absence of an ex- 
ternal magnetic field, the support of a non-trivial P(q s ) 
should range from q s = to <7ea- In Fig. [15] we plot 
P{q s ) at D = and size L = 16 for all simulated tem- 
peratures: as T decreases P(q s ) moves from a Gaussian 
to a bimodal distribution. The important issue is, then, 
whether the continuous part in between the two peaks 
at low T goes to zero or not as L increases. In Fig. [TBI 
we plot P(q s ) at the lowest thcrmalized temperature for 
L = 6,8, 10, 12 and 16 and, in the inset, we plot the 
values of -Pl(O) displaying no decreasing trend with in- 
creasing L. The states, thus, appear to be many and 
different among themselves, since they are found with a 
finite probability within a non-zero continuous range of 
overlap values, including q s = 0. 

Also P(qi) appears to develop a second peak at small 
qi as L increases, and this signature becomes clearer and 
clearer at low temperature as L increases, cf. Fig. Q~7] 
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FIG. 15: Behavior of the overlap distribution P(q) through 
the second order phase transition and in the low-temperature 
phase for L = 16. 




FIG. 17: Link overlap distribution P L (qi) at T = 0.5, D = 
for L = 6,8, 10, 12, 16. Inset: Variance vs. 1/L tends to a 
very small value a 2 qi — 0.0010(7) as L — > oc interpolating 
with a power-law (1.5(1)). 
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implies the non-triviality of the latter. This can be re- 
alized by recalling that in the SK model qi = q 2 and by 
comparing P(qi) to the distribution Q(q a ) of an auxiliary 
variable 



A + Bq 2 s + z^/T~, 



(29) 



with z a Gaussian random variable of variance a z and 
zero mean, that mimics the presence of fluctuations due 
to the finite size of the considered systems. 

At a given point of the phase diagram D,T and for 
a given size L, the parameters A(L), B{L) and a z (L) 
can be obtained by minimizing the Kullback-Leibler 
divergence^ (KLD) between P(qi) and Q{q a )- 



N h 



Dkl[P,Q] = ]TPfe) log 



(30) 



The analysis of FSS of the variance of P{qi) might help 
to better evaluate the breadth of the distribution in the 
thermodynamic limit. Its behavior for various sizes is 
exemplified in the inset of Fig. [17] at the lowest T/T c 
we simulated for D = 0. The variance tends to a small 
finite value and we cannot make a definitive statement 
about P(qi) tending towards a delta distribution, as con- 
jectured by the TNT scenario. Moreover, the study of the 
variance does not yield any indication about the shape of 
the distribution. In particular, about the FSS behavior 
of the two peaks expected in RSB theory. 



Equivalence of site and link overlap 
distributions 



We can, then, implement a more refined analysis of the 
pdf data and check whether P(q s ) and P(qi) are actually 
equivalent and, thus, if the non-triviality of the former 



We will refer to this one as the "left" KLD. The "right" 
KLD is the same formula exchanging P and Q, where 
the symmetrized divergence (sKLD) between P{qi) and 
Q{q a ) is defined as^ 



D Kh [P,Q] 



N bin 

E 

i=l 



Pfe)log 



Qfe) 



Q(*)iog 



Q(Qi) 



P(Qi). 
(31) 

In Fig. [18] we plot, the finite size values of the param- 
eters A and B. Besides the values of the parameters 
minimizing the symmetrized KLD, Eq. (|31[) we also plot 
the values of A and B minimizing the left and the right 
unsymmetrized KLD's. We observe that, as L increases 
the spread between different estimates tends to vanish. 
The infinite size limit of a z is always compatible with 
zero, signaling that FS effects actually tend to vanish as 
L increases, though with large statistical errors at low 
temperature, implying that smaller sizes might hinder a 
correct FSS. 
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As instances we plot the matching of the two distribu- 
tions Q(q a ) and P(qi) in Figs. [19] at T = 0.5 ~ 0.5T C and 
T = 0.7 ~ 0.7T C at size L = 16 and D = 0. In the in- 
sets we plot the size behavior of A and B from the sKLD 
for the two specific cases. In the first case, performing a 
power-law FSS scaling to L — > oo we obtain that B inter- 
polates a negative value! In the second case the L — > oo 
limit yields a positive value. This observation is contrast- 
ing from the behavior, cf. bottom panel of Fig. [TBI of 
B(T) growing with decreasing T at all fixed sizes. Quite 
evidently, the low L strong fluctuations strongly bias the 
interpolation at small T. To show it in a clearer way, in 
Fig. [20] we plot the asymptotic values of both A and B 
for all simulated temperatures both from the sKLD and 
as the average of the extrapolation of the values mini- 
mizing the right and left unsymmetrized KLD's. With 
Aca{T) the two estimates appear to be consistent at all 
temperature and reproduce the qualitative behavior de- 
tected in all finite L cases, compare with Fig. [THJ For 
Boo(T), at low T the two estimates are not consistent 
anymore. Moreover, BooiT) decreases with T below a 
certain T ~ 0.7, unlike its finite L counterparts (at least 
as L > 10), cf. Fig. HFJ 

We face strong finite size effects and a crossover be- 
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tween small and large sizes is taking place. However, due 
to the fact that we cannot easily thermalize larger sys- 
tems at low temperature, we cannot make any definite 
statement on the behavior of Boo (T) for very low T. We 
simply do not have enough reliable points in L at our 
disposal. The finite size behaviors, though, strongly sug- 
gest that Q(q a ) and P{qi) are, indeed, equivalent even 
below T — 0.7. In any case, the equivalence is proven 
for T > 0.7 implying that not only the equilibrium states 
have a non-trivial distribution but also their excitations, 
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yielding evidence in favor of the third scenario consid- 
ered, the RSB theory, rigorously valid in mean-field sys- 
tems. 



B. Position Space Four Spins Correlations 

We now investigate the behavior of the four spins cor- 
relation function, defined in Eq. (fTU|). in position space. 
We recall that the droplet and TNT theories predict that 
Ci{x) tends to a plateau of height q^ A (cf. Sec. IVII[) 
whereas RSB theory predicts for C^x) at T < T c a 
power-law decay ~ x~ a . We, thus, have to compare our 
data with the prediction of one of these hypotheses. 

Since we are dealing with small systems, we must first 
consider possible FS effects. Indeed, because of the peri- 
odic boundary conditions imposed on the simulated sys- 
tem, the correlation function that we actually measure 
at a distance x also contains the contribution of correla- 
tions at distance x + kL, with k — 1, . . . , oo and the true 
(yet unknown) correlation function Ci{x,y,z) is related 
to the measured one - C^{x, y, z) - by the relationship: 



d(x,y,z) = 



0,00 



k-r-L, 



k y L,z + k z L) (32) 



For large distances, when C4 is smaller, these extra con- 
tribution will strongly bias the estimate of the true C4 
behavior in space. In particular, correlations at larger 
distances, of order L/2, will experience relatively stronger 
systematic errors than C4(|r| <C L). 

We will now present our results for the case D = 0. 
For temperatures down to the critical region we simu- 
lated lattices with sides of length up to L = 24. The 
largest thermalized size for T down to 0.5T C is, instead, 
L = 16. In Fig. [51] we plot the x behavior at T = 1.5 
in a log-log plot for the sizes 10, 12, 16, 20, 24. One can 
observe that FS effects are limited to the last point at 
L/2. The rest of the curves completely superimpose. 

At high temperature, correlations are expected to de- 
cay exponentially at large enough distances. As tem- 
perature is lowered towards criticality the C^{x) should 
become power-law eventually decaying as x ~ d+2 ~ v at 
T = T c . We, then, interpolate the four-spins correlation 
function along the x-axis at criticality with the function: 



Cf{x) 



-So 



Jx/i 



-1/5 



(33) 



and equivalently for y and z, due to the anisotropy of 
the system in absence of an external field. This is a 
function containing a crossover between a short distance 
power-law decay, x~ a , and an exponential decay, with 
characteristic 'correlation' length £. In Fig. [5T]the func- 
tion interpolating the L = 24 C±{x, 0,0) is plotted with 
a = 0.402(9), 5 = 0.69(1)/ = 1.25(1) with X 2 = 0.088. 
As the temperature decreases the correlation length in- 
creases until it becomes too long to be observed in the 
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FIG. 21: Correlation between local overlap for D = and 
sizes L — 10, 12, 16, 20, 24, at the largest simulated tempera- 
ture T = 1.5. The fit with Eq. <(33j is also plotted. 
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FIG. 22: Behavior of C 4 {x) at D = 0, for L = 10, 12, 16, 20, 24 
and T — 1 ~ T c . The interpolation with a simple power- 
law, a = 0.64(1), is shown for L — 24. On shorter systems: 
a = 0.64(1), L = 20 and ol = 0.65(2), L = 16. 



analyzed systems. In the inset of Fig. [21] we plot the T 
behavior of t, a and 5 until the fit becomes inconsistent 
T ~ 1.15. 

In Fig. [22] we plot the C4 curves at T ~ T c for sizes 
L = 10, 12, 16, 20, 24, as well as the interpolation of the 
latter with Ax~ a (the correlation length is too long to de- 
tect the exponential contribution in Eq. (|33[) ) . The expo- 
nent equals the power at criticality a = d—2+rj = 0.64(1) 
(at crystal field D = it was 77 = 0.36(1), cf. Tab. HIT]) . 
At T = 1 the interpolated value of a for the L = 24 
Ci(x) curve is a = 0.64(1), a = 0.65(2) for L = 16 and 
a = 0.64(1) for L = 20. FS effects appear to be stronger 
now w.r.t. Fig. [2T1 and evident also for x < L/2 (only 
points for x < L/4 actually stay on the x~ a curve). 

Approaching T c , as T < 1.2, cf. Fig. [23] it is not 
possible to detect a crossover between power-law and ex- 
ponential decay and the simple power-law decay is tested. 
In the inset the power behavior in T is shown and com- 
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FIG. 23: Behavior of C 4 (x) for L — 24 and T = 
0.95, 1, 1.05, 1, 1.15. The interpolation with a simple power- 
law is also shown for L — 24. Inset: behavior of the power 
a vs. T. The dashed vertical and horizontal lines denote, re- 
spectively the estimates of —d + 2 — rj and T c (with errors, cf. 
Tab. HV] T c = 1.01(1), n = -0.35(1)). 
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FIG. 24: Behavior of C±{x) for L = 16 and T = 
0.5,0.6,0.7,0.8,0.9.1.0. The interpolation with a simple 
power-law is also shown. On the right hand side the values of 
Qea(T) are displayed. Inset: behavior of the power a vs. T 
(full line) compared to the a(T) behavior around criticality 
for L = 24 (dotted line). 



pared with the power at criticality, a = 0.65(1). 

Decreasing further the temperature we show in Fig. [24] 
that the behavior is power-law until x ~ L/4 is reached. 
At that point the curves bend upwards as it did at crit- 
icality and even at high temperature, cf. Fig. 1211 This 
bending is, however, an artifact due to the contributions 
induced by the periodic boundary conditions. In Fig. 1241 
on the right hand side, we show the values of q^ A at the 
same temperatures of the plotted C±{x). At all temper- 
atures the C/i{x) soon decays below the corresponding 
value of g| A . For the sizes simulated our data are, thus, 
not consistent with the observation of a plateau at q^ A 
as predicted by the droplet and TNT theories. 



VIII. CONCLUSIONS 

In the present work we have performed Parallel Tem- 
pering Monte Carlo numerical simulations in the temper- 
ature/crystal field plane of the random 3D Blume-Capel 
model on a cubic lattice. This is a spin-1 spin glass, 
whose constituent features try to capture at least one 
supposed mechanism underlying inverse transitions: the 
raise of inactive components at low T. 

In particular, we have analyzed the second order phase 
transition carrying out the computation of the critical 
temperatures and indeces by means of parallel temper- 
ing simulations in temperature at different values of the 
chemical potential D. In this analysis we have carefully 
checked FS effects, identified eventual crossovers from 
small to large size scaling and neglected data for corre- 
spondingly too small sizes. We verified that for different 
values of D the system is always in the same universality 
class (as far as a continuous transition occurs) looking, 
e.g., at different universal scaling functions of £ C /L, such 
as the Binder parameter g, or the quotients of xsg, £c 
and g between systems at L and 2L. The outcome is that 
at all D < D^c the second order transition belongs to 
the same universality class of the 3D Edwards- Anderson 
model for spin-glasses. 

We, then, estimated the position of a tricritical point, 
D ~ 2.1, T ~ 0.5, beyond which the transition is first 
order with jump in density and in overlap parameters. 
This transition is first order in the thermodynamic sense, 
i.e., latent heat is exchanged and, even though the sys- 
tem is disordered, it is not related to the random first 
order transition taking place in structural glasses^ We 
employed and compared four different methods to infer 
the critical line from FS data. This observation con- 
firms the claim of Fernandez et al. 56 about the existence 
of such transitions in quenched disordered short-range 
finite-dimensional systems. In the present model the first 
order transition can be seen by means of standard paral- 
lel tempering algorithm in the canonical ensemble, simply 
tuning an external pressure-like parameter. 

The first order transition line has the property of dis- 
playing inverse freezing, as can be observed from the 
phase diagram, cf. Figs. HH [14] the low temperature 
phase is paramagnetic and the system 'freezes' into a 
spin-glass phase as T is increased. This is at difference 
with the thermodynamic behavior of the original, or- 
dered, BC model (mean- field or finite dimensional) j 44 i 51 
In presence of quenched disorder, a low temperature 
paramagnetic phase exists that can acquire a very low 
density and this is the source of the entropy decrease with 
respect to the high temperature paramagnetic phase. 

Both the inverse freezing transition and its first order 
nature were not observed in the same model on a hierar- 
chical latticed 

Eventually we present our analysis of the overlap dis- 
tribution functions and the four-spins correlation func- 
tions at criticality and in the glassy phase, at D = 
for T down to 0.5T C . From the behavior of site overlap 
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distribution at zero overlap, Pl{<1s — 0), and from the 
variance of the link-overlap distribution Pl(Qi) we get 
evidence in favor of a complex organization of states in 
the SG phase, displaying features typical of the Replica 
Symmetry Breaking theory holding for mean-field sys- 
tems (d > 6). We cross-checked this observation com- 
paring, with the Kullback-Leibler divergence, the link- 
overlap distribution with the distribution of a function 
of the squared site overlap, q a ~ A + Bq 2 s . We carefully 
analyzed the finite size effects at low temperature finding 
that for T < 0.7 small size fluctuations strongly bias our 
estimates, yielding negative B coefficients of the term, 
decreasing with temperature, unlike any finite size B{T) 
behavior. In order to have a self-consistent estimate we 
would need to thermalize at T > Q.7T C systems of size 
sensitively larger than L = 16. 

Looking at the position dependence of the four-spins 
correlation functions we are able to detect, for T > 1.2T C , 
a crossover between a short-distance power-law decay and 
a long-distance exponential decay and we can identify a 
length-like parameter £ playing the role of the correlation 
distance, growing as T decreases. As the critical temper- 
ature is approached and I becomes similar to the maxi- 
mum feasible distance in the simulated system (~ L/2), 



Ci{x) can be interpolated with a simple power-law. We 
checked that for sizes L — 16, 20 and 24 the exponent of 
C±(x) at T c is equal to d — 2 + rj, where n = —0.36(1) 
is the value obtained from the analysis of the critical 
properties performed with the quotient method. We also 
probed the power-law behavior for temperatures down 
to 0.5T C at distances far away from border, where finite 
size correction are too strong. Indeed, periodic boundary 
conditions systematically increase correlations, above all 
where they are small (or vanishing), i.e., at large dis- 
tance. We compare the low temperature behavior with 
the prediction of TNT and droplet theories that C±(x) 
should tend to a plateau C4 ~ q| A for large x. Even 
though we are not able to reach "large x n , we show that 
Ci{x) < gj| A already at small distance. 
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